In [165]:
import numpy
import matplotlib
import numpy as np
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
%matplotlib inline
from pylab import *
Liddriven Cavity :
Before exploring this example note that we changed the basic code and removed loops . This really helps us to save the time . If you run this code with considering loop instead of indexes it will take approximately 15 hours to run .
In this simulation we have a 2D fluid flow that is driven by a lid at the top which moves at a speed of uo . the cavity has four sides which three of them are solid and no-slip boundary condition is considered for walls and the top one which has moving velocity . If you remember in the file (illustration) we mention how to deal with moving or solid boundaries . This problem is twodimensional .
We need to consider both advection and diffusion roles in Navier-stokes to ideally simulate this problem . In this example in contrast to previous examples the fluid has the ability to move .
Question : How can we measure the strength of motion or strength of the fluid induced by moving top lid ? We have the answer thanks to the fluid mechanic course . The fluid could be Water or Gas with different kinematic velocities or densities .
We at first need to non-dimensionalize the navier stokes equation . If we do the Reynolds number will appear which is defined as the inertia force over viscosity force . Higher reynolds number implies higher inertia or lower viscosity resistance and vice verca .
What changes should be considered to shift our code from conduction to fully Navier-Stokes equation ?
1- define parameter reynolds which you want to consider , It means that you definitely understand that you can not directly the parameters which you directly use in FDM or in original format of your equation . For example if you set the reynolds number with adjusting the top velocity or kinematic viscosity in original equation here you are not allowed to consider any arbitary value .
The most important point is that you must make equal value of Reynold number . You need to have reynold number for example 1000 then you choose the charectristic velocity as 0.1 ( large value breaks low mach number ) and also you need to have relaxation time lower than 0.5 ( Chapman-Enskog for navier stokes suggestes us ) . We extract the viscosity from reynolds definition after adjusting characteristic velocity ( lid velocity) and characteristic length (ny- the number of lattice nodes in y direction) . By obtaning kinematic viscosity relaxation time could be obtained by satisfying the mentioned criteria .
2- We already set D2Q9 lattice velocities so we do not need to change streaming section however collision distribution function should be modified to account the advection roles . ( Make attention how we change this par)
ezafe kardane formule vase ghesmate collision
3- comparison were made by "Ghia.et al" . we used the u velocity along the horizontal line and V velocity along vertical line at the middle of the cavity .We got these data from " Gerris Flow Solver " site .
Further efforts :
1- Try to change the Reynolds number to see what would be the vortex behavior . 2- Try to change the number of nodes in Lattice . what is the relation between lattice nodes and the length of the cavity? 3- Try to change the lid driven position ( for example change the left side from no-slip to moving boundary) 4- Try to creat cavity with unequal sides' length . Lets the length of the cavity is twice of the width .
In [166]:
nx=100 # the number of nodes in x direction lattice direction
ny=100 # the number of nodes in y direction lattice direction
Re=1000 # rayleigh number
uo=0.1
rhoo=1.0
visco=uo*ny/Re
D=3.0 # the dimension of the problem
omega=1./(0.5+(visco*D)) #Chapman-Enskog dt=1 and dx=1 relaxation flow
mstep=100000 # the number of time step
k=9 # k=0,1,2,3,4,5,6,8,9
In [167]:
w=numpy.ones(k) # witghting factor
cx=numpy.ones(k)
cy=numpy.ones(k)
rho=numpy.ones((ny+1,nx+1) ) # density matrix
u=numpy.ones((ny+1,nx+1) )
v=numpy.ones((ny+1,nx+1) )
f= numpy.ones((k, ny+1,nx+1)) # distribution function
feq= numpy.ones((k, ny+1,nx+1))
t1=numpy.ones((ny+1,nx+1) )
t2=numpy.ones((ny+1,nx+1) )
rhon=numpy.ones(nx+1)
In [168]:
y1= numpy.zeros(17)
u1= numpy.zeros(17)
y1[0]=-0.49933 ;u1[0]= -0.000882
y1[1]=-0.444335 ;u1[1]= -0.181701
y1[2]=-0.43629 ;u1[2]= -0.201989
y1[3]=-0.428914 ;u1[3]= -0.222276
y1[4]=-0.397406 ;u1[4]= -0.297251
y1[5]=-0.327052 ;u1[5]= -0.383699
y1[6]=-0.217948 ;u1[6]= -0.27788
y1[7]=-0.046595 ;u1[7]= -0.106804
y1[8]=0.001598 ;u1[8]= -0.060949
y1[9]=0.118733 ;u1[9]= 0.057217
y1[10]=0.235193 ;u1[10]= 0.186849
y1[11]=0.352315 ;u1[11]= 0.333239
y1[12]=0.45404 ;u1[12]= 0.466401
y1[13]=0.461386 ;u1[13]= 0.511382
y1[14]=0.469392 ;u1[14]= 0.574884
y1[15]=0.476719 ;u1[15]= 0.659554
y1[16]=0.5 ;u1[16]= 0.999118
In [169]:
x1= numpy.zeros(17)
v1= numpy.zeros(17)
x1[0]=-0.5 ;v1[0]=0.00069404
x1[1]=-0.43768 ;v1[1]=0.275621
x1[2]=-0.429602 ;v1[2]=0.290847
x1[3]=-0.421523 ;v1[3]=0.303994
x1[4]=-0.406521 ;v1[4]=0.326826
x1[5]=-0.343624 ;v1[5]=0.371038
x1[6]=-0.273803 ;v1[6]=0.330015
x1[7]=-0.265724 ;v1[7]=0.32307
x1[8]=-0.000289 ;v1[8]=0.0252893
x1[9]=0.304962 ;v1[9]=-0.318994
x1[10]=0.359781 ;v1[10]=-0.427191
x1[11]=0.40652 ;v1[11]=-0.515279
x1[12]=0.445182 ;v1[12]=-0.392034
x1[13]=0.45326 ;v1[13]=-0.336623
x1[14]=0.461339 ;v1[14]=-0.277749
x1[15]=0.46884 ;v1[15]=-0.214023
x1[16]=0.5 ;v1[16]=-6.20706e-17
In [170]:
strf=numpy.ones((ny+1,nx+1) ) # streamfunction
##================ Initial boundary condition
w[0]=4./9 #4./9
w[1:5]=1./9.
w[5:9]=1./36.
cx[0]=0.0
cx[1]=1.0
cx[2]=0.0
cx[3]=-1.0
cx[4]=0.0
cx[5]=1.0
cx[6]=-1.0
cx[7]=-1.0
cx[8]=1.0
cy[0]=0.0
cy[1]=0.0
cy[2]=1.0
cy[3]=0.0
cy[4]=-1.0
cy[5]=1.0
cy[6]=1.0
cy[7]=-1.0
cy[8]=-1.0
In [171]:
##================== Initial value
rho[0:ny+1,0:nx+1]=rhoo
v[0:ny+1,0:nx+1]=0.0
u[0:ny+1,0:nx+1]=0.0
u[ny,0:nx+1]=uo
for i in range(nx+1):
for j in range(ny+1):
for l in range (k): #k=0,1,2,3,4
f[l,j,i]=w[l]*rho[j,i]
In [172]:
## Main loop : comprised two parts :collision and streaming
##=====================
for n in range(mstep) :
#print n ,u[nx/2,ny/2]
if (n%10000==0):
print n
# collision process fluid flow
t1[0:ny+1,0:nx+1]=u[0:ny+1,0:nx+1]*u[0:ny+1,0:nx+1]+v[0:ny+1,0:nx+1]*v[0:ny+1,0:nx+1]
for l in range (k):
t2[0:ny+1,0:nx+1]=cx[l]*u[0:ny+1,0:nx+1]+cy[l]*v[0:ny+1,0:nx+1]
feq[l,0:ny+1,0:nx+1]=w[l]*rho[0:ny+1,0:nx+1]*( 1.+ 3.*t2[0:ny+1,0:nx+1]+4.5*t2[0:ny+1,0:nx+1]**2-1.5*t1[0:ny+1,0:nx+1] )
f[l,0:ny+1,0:nx+1]=(1.-omega)*f[l,0:ny+1,0:nx+1]+omega*feq[l,0:ny+1,0:nx+1]
f[1,0:ny+1,nx:0:-1]=f[1,0:ny+1,nx-1::-1]
f[3,0:ny+1,0:nx]=f[3,0:ny+1,1:nx+1]
f[2,ny:0:-1,0:nx+1]=f[2,ny-1::-1,0:nx+1]
f[5,ny:0:-1,nx:0:-1]=f[5,ny-1::-1,nx-1::-1]
f[6,ny:0:-1,0:nx]=f[6,ny-1::-1,1:nx+1]
f[4,0:ny,0:nx+1]=f[4,1:ny+1,0:nx+1]
f[8,0:ny,nx:0:-1]=f[8,1:ny+1,nx-1::-1]
f[7,0:ny,0:nx]=f[7,1:ny+1,1:nx+1]
## =============================
#left Boundary- Bounce back
f[1,0:ny+1,0]=f[3,0:ny+1,0]
f[5,0:ny+1,0]=f[7,0:ny+1,0]
f[8,0:ny+1,0]=f[6,0:ny+1,0]
#right Boundary-bounce back
f[3,0:ny+1,nx]=f[1,0:ny+1,nx]
f[7,0:ny+1,nx]=f[5,0:ny+1,nx]
f[6,0:ny+1,nx]=f[8,0:ny+1,nx]
# top Boundary -bounce back
rhon[0:nx+1]=f[0,ny,0:nx+1]+f[1,ny,0:nx+1]+f[3,ny,0:nx+1]+2*(f[2,ny,0:nx+1]+f[6,ny,0:nx+1]+f[5,ny,0:nx+1])
f[4,ny,0:nx+1]=f[2,ny,0:nx+1]
f[7,ny,0:nx+1]=f[5,ny,0:nx+1]-(rhon[0:nx+1]*uo/6.0)
f[8,ny,0:nx+1]=f[6,ny,0:nx+1]+(rhon[0:nx+1]*uo/6.0)
# bottom Boundary -bunce back
f[2,0,0:nx+1]=f[4,0,0:nx+1]
f[5,0,0:nx+1]=f[7,0,0:nx+1]
f[6,0,0:nx+1]=f[8,0,0:nx+1]
for i in range(nx+1):
for j in range(ny+1):
sumr=0.0
for l in range (k):
sumr=sumr+f[l,j,i]
rho[j,i]=sumr
for i in range(1,nx):
for j in range(1,ny):
usum=0.0
vsum=0.0
for l in range (k):
usum=usum+f[l,j,i]*cx[l]
vsum=vsum+f[l,j,i]*cy[l]
u[j,i]=usum/rho[j,i]
v[j,i]=vsum/rho[j,i]
u[:,0]=0.
u[:,nx]=0.
u[0,:]=0.
u[ny,:]=uo
v[:,0]=0.
v[:,nx]=0.
v[0,:]=0.
v[ny,:]=0.
In [173]:
strf[0,0]=0.0
for i in range(nx+1):
rhoav=0.5*(rho[0,i-1]+rho[0,i])
if(i>0) :
strf[0,i]=strf[0,i-1]-rhoav*0.5*(v[0,i-1]+v[0,i] )
for j in range(1,ny+1):
rhom=0.5*(rho[j,i]+rho[j-1,i])
strf[j,i]=strf[j-1,i]+rhom*0.5*(u[j,i]+u[j-1,i])
x=numpy.linspace(0,nx,nx+1)
y=numpy.linspace(0,ny,ny+1)
In [174]:
mx, my = numpy.meshgrid(x,y)
plt.figure(figsize=(8,5))
plt.contourf(x,y,strf,20)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.colorbar();
In [190]:
mx, my = numpy.meshgrid(x,y)
plt.figure(figsize=(8,5))
plt.contour(x,y,strf,900)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.colorbar();
In [175]:
mx, my = numpy.meshgrid(x,y)
plt.figure(figsize=(8,5))
plt.contourf(x,y,u,20)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.colorbar();
In [176]:
mx, my = numpy.meshgrid(x,y)
plt.figure(figsize=(8,5))
plt.contourf(x,y,v,20)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.colorbar();
In [177]:
x1=(x1*ny)+(ny/2)
y1=(y1*ny)+(ny/2)
vv=v[ny/2,0:nx+1]/uo
uu=u[0:ny+1,nx/2]/uo
plt.plot(x1,v1, 'r*',label=' Ghia et.al');
plt.plot(x,vv, 'b-',label=' LBM ');
plt.legend();
In [178]:
plt.plot(y1,u1, 'r*',label=' Ghia et.al');
plt.plot(y,uu, 'b-',label=' LBM ');
plt.legend();